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ABSTRACT 

A technique to study combined influence of environmental and genetic factors on the base 
of changes in phenotype distributions is presented. Histograms are exploited as base 
analyzed characteristics. A continuous time, discrete state Markov process with piece- 
wise constant interstate transition rates is associated with evolution of each histogram. 

The technique was applied to IQ longitudinal data (6- and 14-year children) and made it 
possible to draw conclusions concerning development of Russian children before school 
and during the first stage of school education as well as dependence of environmental and 
genetic influence on IQ level. 


INTRODUCTION 

Combination of genetic and environment factors is a logical consequence of the 
grouping humans into population of self-determining individuals who share both genes and 
environment in common. Combined genotype-environment effects conditioned by assortative 
mating, genotype-environment correlation, genotype-environment interaction, etc. are 
referred to here, as a whole, as the results of influence of some systematic factor. 

Longitudinal analysis of such influences is frequently of interest in psychological 
investigations. 

Means for this analysis are usually include various tests for concordance, factor 
analysis of variance/covariance structures and, as one of its most advance variants, analysis of 
simplex models (Guttman, 1954; Boosma, Martin & Molenaar, 1989; Dolan, Molenaar & 
Boosma, 1991). Nevertheless, analyzing such integrating characteristic as variances, 
covariances and means, a researcher loses a lot of information. As a rule, he also has 
difficulties during extrapolation of obtained results to the time domains between observation 
points and during study of dependences of factor effects on different levels of individual 
characteristics. That is why the search of other solution techniques based on new principles is 
topical. 

Presented here is the technique intended for studying evolution of measured 
phenotype distributions with age. It follows the ideas stated in our previous paper (Kuravsky 
& Malykh, 2000) where genetic differences between twins and effects of non-shared 
environment were estimated. As earlier, histograms are used as base analyzed characteristics. 
However now they describe individual characteristic distributions for unrelated individuals, 
but not the differences in twin pairs. Using measurements in a series of time points, the 
technique discussed makes it possible to determine how the systematic factor changes an 
initial phenotype distribution during the selected time period. 
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The technique is demonstrated on the analysis of IQ longitudinal data (6- and 14-year 
children). The main purpose of this work was to investigate combined environmental and 
genetic influences on population as functions of age and IQ level. Some interesting 
conclusions concerning development of Russian children before school and during the first 
stage of education were drawn from the analysis. In particular, it was shown that standard 
education promoted development of unable individuals more effectively than development of 
capable ones. Results of analysis on the base of the presented approach were compared with 
the output obtained with the aid of classical phenotypic simplex model. Finally, a so-called 2- 
dimensional model to compare combined genotype-environment effects leading to changes in 
population phenotype and effects of non-shared environment evaluated on twin data is 
discussed. 

Previous research covering similar areas include Markov chain Monte Carlo methods 
(Eaves & Erkanli, 2003; Creutz, 1979; Gelfand & Smith, 1990) and random utility models 
(Falmagne & Regenwetter, 1996). Both approaches carry out histogram estimation. The first 
one constructs a Markov chain on the parameter space of unknown quantities such that, 
starting with a series of trial values (e. g. means, regressions, etc.), after an initial series of 
iterations, successive iterations represent samples from the unknown joint stationary 
distribution. The Gibbs sampler (Creutz, 1979; Gelfand & Smith, 1990) is popular approach 
of this type to construct a Markov chain with desired properties. As distinct from the 
approach presented in this paper, Markov chain Monte Carlo methods are not acceptable for 
longitudinal analysis of observed data because of both their inefficient representation of 
evolutionary processes and inability to solve the inverse problem, viz.: to derive necessary 
development characteristics from the observed distributions. 

The second approach is represented by a model for approval voting. It is based on the 
notion that each voter has a personal ranking of the alternatives, with a random utilities being 
defined by the probability distribution on the set of all rankings. However, this method is 
overly attached to the details of a particular application. Its model is static and, in contrast to 
the technique presented in this paper, is not acceptable for representation of evolutionary 
processes. Moreover, the random utility model fails in case of violations of the basic axiom 
system (the authors showed a corresponding example). 

METHOD 

In the model presented, the distribution of individual characteristics at the moment of 
birth describes influences of heritability and perinatal/prenatal environmental factors. As a 
result of subsequent influences of the systematic factor, initial distributions are transformed to 
current ones during the time of observation. This idealization is based on the results of twin 
studies of IQ, which demonstrated that both environment and genetic influences on individual 
differences are present and dynamic throughout development (Eaves, Eysenck & Martin, 

1989; Malykh, Egorova & Meshkova, 1998). As genetically uninformative samples are 
analyzed, effects of these factors are not separated here. 

It is important for the parameter of interest to be measured in a time- independent 
scale (i. e. the testing technique must take into account the natural age development as is the 
case for IQ). 

A continuous time discrete state Markov process with piece-wise constant interstate 
transition rates 1 is associated with evolution of each histogram. Given a continuous trait, the 
available actual range of the analyzed quantity is divided into several intervals. These 
intervals are considered as separate discrete states in which the trait value has some 
probability to find itself. In due course transitions between the states are the case. The number 
of these states is determined by desirable precision of estimates and available sample size 2 . 


-j 

These processes are frequently used to solve problems of the queuing theory. 

2 Loss of information due to discretization is not a serious problem, especially if large enough samples 
are possible such that state intervals are smaller than measurement error. 
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Some considerations concerning variation of state interval length are presented in our 
previous paper (Kuravsky & Malykh, 2003). 

This model may be represented by a graph (Figure 1) in which nodes (depicted as 
rectangles) correspond to the states; branches (depicted as arrows) correspond to transitions. 
The trait evolution process may be imagined as a random walk along the graph from one state 
to another following the arrows. Time is supposed to be continuous. State-to-state transitions 
are instantaneous and take place at random time points. These transitions are effects of the 
systematic factor. To describe it mathematically means to show how this factor has an 
influence on population. 



/// Pj jUj+1 Hn 

Figure 1. The model used to describe evolution of phenotype distributions with age, where 
the Xj ( j = 0,1,..., /I ) are states, the Xp (p — 0,1,. ..,w 1) and Mq 

( q = 1 , 2 ,..., n ) are flow rates. 


It is assumed that state-to-state transitions (corresponding to each branch of the 
graph) meet the properties of Poisson’s flows of events (Ovcharov, 1969). It may be proved 
(Ovcharov, 1969) that the number of events X in these flows falling into any interval of the 
length T adjoining to time point t is described by the Poisson distribution: 


P t ^(X = m) = 


ml 


in 


e 5 


where P t T (X = m) is the probability of appearance of 111 events during the considered 
interval, a(t,r) - mean number of events falling into an interval of the length T adjoining 

to time point t . Only stationary flows (where a(t,r) = Tjr , 77 = const ) will be taken 
up here. Parameter TJ is the rate of a stationary flow. It is equal to mean number of events per 
unit time interval. Mean time interval between two adjacent events in this flow is \l 1J . 

The aforementioned assumptions concerning the nature of event flows are usual for 
applications as these flows (or flows which close to them) are frequently take place in reality 
because of limit theorems for events (Ovcharov, 1969). 

The system shown in Figure 1 is a finite chain of n + 1 states where transitions from 
the state X j ( j = 1 , 2 ,..., n — 1 ) are possible only to the preceding state Xj_\ or to the 

next state . Available from the states Xq and X n are the states Xj and X n _\ 
correspondingly. 

Flow rates are denoted as Xp and /Uq ( p — 0,1,..., n — 1 and q = 1,2,.. .,71 ). 
Parameters Xp define how the influences of the systematic factor under study promote 

increase of the examined phenotype values. Parameters Pq describe the inverse effect of this 
factor - decrease of the phenotype. 
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If one denotes possible bottom and top limits of parameter change as Dbot and 

n n n Dfop ~ Dbot 

Is top > the state Xq corresponds to the interval from Db ot to Db 0 t d , the 

1 n + 1 

, • , r\ D tOp ~ D bot _ _ Dfop ~ D bo t 

state X] - to the interval from Dbot d to D bot ^ , and 

n + 1 n + 1 

so on. The following set of ordinary differential equations (Ocharov, 1969) may be drawn to 
describe the time history of state probabilities: 

dpo (t ) 


dt 


-AoPoiO + Mi p x ( 0 ; 


dp jit) 


dt 


(Aj + Mj)Pj(t) + Aj_ x pj_ x {t) + jU J+l p J+l (t); 


= -PnPniO + K-\Pn-\it) 

dt 

where Pk ( t ) is the probability to be within the state Xj ( at the time point t 

(k = 0,1,..., n ). 

To integrate these equations, one has to assign initial conditions 

n 

Po (0), p\ (0),..., p n (0); 2>*(°) = i- 

k—0 

n 

The normalization condition Y;Pk( t ) = 1 is valid at any time point. It is postulated that 

A=0 

t = 0 is the moment of birth. 

Since the quantities that are formed as a result of influence of many different 
elementary factors (such as IQ) may be considered as asymptotically normal ones, it is 
supposed that at the moment of birth and other test time points of interest, these 
characteristics may be approximately described by some normal distributions. Mean HI q and 

standard deviation CTq at point t = 0 characterize the initial distribution. 

So, estimation of systematic effects are brought to the calculation of mean HI q , 

standard deviation CTq and flow rates A p and Pq . Values ensuring the best fit of expected 

and observed frequencies of falling into each state at the specified time points are taken as 
estimations of these free parameters. Estimated values are considered as the characteristics of 
systematic effects which have become apparent during observations. Normality of phenotype 
distributions at the moments of observations (as is the case for IQ study) is actually an 
additional constraint for selection of flow rates. 

In case of normally distributed phenotypes, free model parameters are to be 
determined as functions of means and standard deviations of corresponding normal 
distributions to avoid model adjustment to disturbances induced by sampling errors. Such 
smoothing of observed histograms reduces the dependence of final conclusions from these 
errors. 
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Expected state probabilities are calculated by means of integration of the presented 
set of differential equations. Expected frequency to fall at the k -th state equals to PkN, 

where - probability of being in this state, IS - number of cases. Corresponding original 


observed frequencies F/ ( result from measurements obtained during longitudinal study. 

Provided that the given expected frequencies describe observed data 3 , the following Pearson 
statistic is distributed asymptotically according to a chi-square distribution (Cramer, 1946): 


I 

k—0 


( F k - Pk N y 


PkN 


One should regard this sum as a goodness-of-fit measure in the sense that its large values 
correspond to bad fit and its small values correspond to good fit. Under certain general 
conditions (Cramer, 1946; Fisher, 1924), the number of degrees of freedom is equal to 
n — l , where / is the total number of free parameters determined from the sample under 
study, and serves as a standard by which one can judge whether such measure is large or 
small. 


Estimations of free parameters are found as the values minimizing the sum of 
goodness-of-fit measures at specified time points in which observed data are available, with 
the original random observed frequencies in the expressions of Pearson statistic being 
replaced by the corresponding frequencies of the best- fitted normal distributions. Such 
substitutions disable the opportunity to use a chi-square distribution in the test for 
concordance. Elowever, one can use sum of the aforementioned statistics simply as a criterion 
for minimization to find the best- fitted values 4 and, then, carry out the test for concordance of 
the obtained model and original observed data. 

The employed estimation procedure consists of three stages. On the first stage, some 
numerical integration scheme for the aforementioned differential equations is coded to 

calculate all Pk it) using the Microsoft® Excel spreadsheet. The probability functions are 

computed with some specified time step h from initial zero time point to the given specified 
upper time bound. Runge-Kutta methods (Bakhvalov, 1975) (or their equivalents) proved to 
be sufficient to get acceptable accuracy of solution. It is of vital importance that Excel 
supports dynamic li nks between cell contents, viz.: if one locates current values of free 
parameters and time step h in the separate cells to which the cells containing formulas for 
calculation of Pk{t) and initial state probabilities are referred, all the solution will be 

automatically modified when values of free parameters are changed. 

On the second stage, parameters of best-fitted normal distributions for observed 
frequencies F \ are estimated, and, later, the frequencies defined by these distributions are 


used instead of the observed ones. Parameters of the best-fitted distributions may be 
calculated either as sample point estimates or by the chi-square minimum approach directly. 

On the third stage, a numerical optimization procedure to get required values of free 
parameters is run. The authors used a macros realizing the procedure of non-linear 
optimization called Generalized Reduced Gradient (GRG2). 

It is necessary to note that this technique based on capabilities of spreadsheets, in 
fact, determines coefficients of model differential equations using a given observed solution 
(the inverse problem is solved). 


3 Null hypotheses. 

4 In fact, we apply the least-squares method instead of the method of chi-square minimum in its 
classical form. The expression of Pearson statistic is used formally as normalized least squares 
criterion. 
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RESULTS 

In the following sections, application of the presented approach to analysis of IQ 
longitudinal data is discussed. Test data for 94 six-year old children and 70 fourteen-year old 
children are analyzed. All the children entered standard Russian schools. Measurements of 
general IQ (GIQ) were carried out twice: at age 6 and age 14 (for the same children). 

Expected frequencies were adjusted to observed ones at the time points corresponding to the 
6- and 14-year age. IQ range from 60 to 160 was taken into account. The model to be fitted 
was represented by a chain containing 20 states. 

Chi-square tests show no significant differences between observed frequencies for 
selected IQ intervals and corresponding expected frequencies of normal distributions (p = 

0.83 for 6-year sample and p = 0.88 for 14-year sample). Therefore the hypothesis of 
normality fits our data. 

According to F-test, there are no significant differences between variances at ages 6 and 
14 (p = 0.12). So we can accept the hypothesis of homogeneity of variances in the following 
considerations and use their common value in estimations. In its turn, /-test shows highly 
significant differences in means between the samples in question (p < 0.00001). 

The sum of goodness-of-fit measures at the 6- and 14-year time points was minimized by 
selecting the following free parameters: 

• flow rates X p and Mq for the range 0-6 years; 

• flowrates X p and JUq for the range 6-14 years; 

• mean HI q and standard deviation <Tq . 

To keep reasonable balance between the model complexity and depth of investigation, we 
have to use minimal number of free parameters that ensures acceptable model fitting to 
observed data and capabilities for analysis of interest. To determine this number, three models 
of different capacity were compared. In the first model, flow rates were assumed to be the 
same for the ranges 60-85, 85-1 10, 1 10-135, 135-160 of IQ units: Xq = X] = X 2 =Xt ) ,..., 

^16 = / ^17 =Xig =X i9 ; H\ — H 2 ~ M3 = M 4 ’ •••> Mil = MlS = Ml9 = / / 20-( Itis 
reasonable to assume that flow rates do not have great changes with GIQ growth.) In the 
second model, these IQ ranges are wider: 60-110 and 110-160. And in the third model flow 
rates do not depend on IQ values. Results of model fitting are given in Table 1. 

Table 1. Fitting models of difference capacity. 


Model 

No 

Goodness-of-fit measure 

(with respect to the 
frequencies of best- fitted 
normal distributions) 

Chi-square goodness 
measure (with respect to the 
original observed data) 

Degrees 

of 

freedom 

P- 

value 

1 

0.43 

19.07 

8 

0.01 

2 

0.79 

18.75 

16 

0.28 

3 

1.75 

18.37 

20 

0.56 


One should select the best-fitted model in the sense of our goodness-of-fit measure, which 
ensures the acceptable correspondence to observed data. Since the first model did not agree 
with them (p=0.01 ), the second model (p=0.28 ) was selected for the following analysis. The 
Third model was rejected as it had the worse value of the minimization criterion. 

In calculations, one year was used as a unit of time measurement. The modified Euler 
method was coded in the spreadsheet as a numerical integration scheme (with time step from 
0.006 to 0.008 years). 
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Figure 2 and Table 2 present obtained estimates of flow rates and IQ initial mean and 
standard deviation. Observed and expected state probabilities for ages 6 and 14 are presented 
in Figure 3. Given in Figure 4 are plots showing the envelopes of distributions of expected 
state probabilities at the moment of birth and test time points ( t = 6, 14). IQ of neonates is 
considered here as an extrapolated characteristic describing their later abilities. It cannot be 
estimated directly. 



Figure 2 . Estimates of flow rates (In - increase; De - decrease). 


Table 2. Model 2: parameters of expected normal distributions. 


Age 

Mean of normal 
distribution 

Standard deviation of 
normal distribution 

0 years 

99.18 

11.81 

6 years 

98.24 

12.63 

14 years 

109.03 

12.63 


To judge the correctness of model application, it is necessary to find out how stable 
are the estimations of free parameters as well as derived conclusions to variations of initial 
data. As to the given problem, such three independent statistics, as sample means for ages 6 
and 14 and the respective common sample standard deviation, together with the important 
assumption on normality of distributions, determine our model. A set of solutions 
corresponding to prescribed initial data confidence intervals was analyzed in paper 1 ' 5 to 
substantiate stability of the obtained results. 
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Figure 3 . Observed and expected state probabilities for ages 6 and 14. 


DISCUSSION 

As a rule, Russian children start their school education at age 6 or 7. Therefore the 
range 0-6 corresponds mainly to preschool development, and the range 7-14 - to the first stage 
of school development. Taking this fact into account, the following qualitative conclusions, 
which are valid within 80% confidence region of the above mentioned independent statistics, 
may be drawn from the results shown in Figure 2 and the stability analysis presented in paper 
(Kuravsky & Malykh, 2003): 

a) during preschool development, the systematic promotion of GIQ decrease is, as a 
rule, greater than the promotion of GIQ increase; 

b) during school development, the systematic promotion of GIQ decrease is, on the 
contrary, much less than the promotion of GIQ increase; 

c) the effect of systematic promotion of GIQ increase during children’s school 
development is, as a rule, much greater than this effect during preschool 
development; 

d) during school development, the GIQ decrease flow rate is negligible for the IQ values 
that are less than the average GIQ; 
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Figure 4. Envelopes of distributions of expected state probabilities at the moment of 
birth and test time points. 


e) during school development, the G1Q decrease flow rate is significant for the IQ 
values that are greater than the average G1Q; 

f) during school development, the G1Q increase flow rate for above the average IQ 
values is less than one for below the average IQs; 

g) during school development, the G1Q decrease flow rate for above the average IQ 
values is much greater than one for below the average IQs; 

h) during preschool development, the increase and decrease scores do not depend 
significantly on the IQ level. 

These conclusions confirm general expectations to the effect that standard education does 
not promote development of capable children. However, the aforementioned /-test (see 
“Results”) shows the efficiency of Russian school education on the average: mean values 
after entering school became significantly greater. At the same time, preschool development 

does not yield essential shift in means: even if the expected value of III q would result from 

direct point estimate by some sample of 1 0 times greater size than the sample size for age 6, 
there would be no statistically significant differences between two means according /-test. 
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Figure 5 . Two-dimensional representation of dynamics of individual characteristics. 
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It is of interest to compare obtained levels of increase and decrease flow rates 
discussed here and levels of convergence and divergence flow rates, which are rode on twin 
non-shared environment and studied in the paper (Gutman, 1954), for the same 5-unit IQ state 
intervals. In fact, it means to consider a 2-dimensional model describing dynamics of two 
independent characteristics (Figure 5). Juxtaposition with data of Table 3 shows that only the 
increased flow rates during school development may be significantly greater than the rates 
induced by twin non-shared environment. Sometimes, decreased and increased flow rates 
during preschool development are comparable with non-shared environment rates. At the 
same time, decrease characteristics during school development are much less. These facts 
demonstrate the important role of random undirected influences on population. 

To appreciate the research effectiveness of the presented method in comparison with 
the allied approaches, the sample under study was examined with the aid of classical 
phenotypic simplex model, which is traditionally used in analysis of repeated measures. 
Following the idea of this model, variables 7ji and t )2 represent observed IQ measurements at 
ages 6 and 14. Their relationship is expressed by the first-order autoregression equation: 

1)2= PlT!l + £>, 

where Pi is the regression of the observed variable at age 14 on the previous variable, 
represents a random input term (innovation 5 ) that is uncorrelated with 1 ) 1 . This simplex model 
is illustrated graphically in the path diagram (Neal & Cardon, 1992; Joreskog, 1970) in Figure 
6 . 



^2 


Figure 6 . Two-point simplex model for studying IQ evolution: variables t)i and 7 ) 2 represent 
observed IQ measurements at ages 6 and 14; Pi is the regression of the observed variable at age 14 
on the previous variable; £ represents an innovation at age 14. 


5 The innovation is that part of the variable at the current time point that is not caused by the variable at 
the previous discrete time point. It differs from random measurement errors that do not influence 
subsequent observed variables. Measurement errors, which are frequently included in the classical 
simplex models, are not considered here. 
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The model in question has three free parameters: rji, £, and /A For two time points of 
observed measurements, the expected covariance matrix of variables rp and 77 ? equals: 

r var 77j (3 2 \dxrj x N 

V A V ar^ Pi var^+var^, • 

Since the number of free parameters is equal to the number of independent statistics in the 
observed covariance matrix, fitting the expected matrix to the observed one yields the perfect 
accordance and the following estimations: var Tji=l 17, var £=722, var /3i=0. 7. 

Both the innovation and regression coefficient are statistically significant: provided 
that each of them taken separately equals to zero, goodness-of-fit measure in the form of 
maximum likelihood function, which is distributed as a chi-square under some general 
conditions (Bollen, 1989), gets statistically significant increase for one degree of freedom. 

Sufficiently great proportion (68%) of observed variance at age 14 is due to 
innovation. The original variance at the first time point explains approximately two times less 
part - 32%. This is all scanty information, which may be obtained from the discussed 
covariance structure. It is much less than we can obtain with the aid of the proposed 
technique. 

Of course, more observation points might be taken and subtler analysis might be 
carried out, but it is well known that acquisition of experimental data is more labor- 
consuming and expensive process than data analysis. Markov model helps us not only to learn 
more about reasons for phenotype changes, but extrapolate observed data to inaccessible time 
points (all points between zero age and age 14). However, using simplex and Markov models, 
we analyze different characteristics: flow rates versus proportions of observed variances 
explained by different factors. Markov models in their current form do not give the 
opportunity to analyze covariance (correlation) structures. 

So, on the same observed data, Markov models, as a rule, yield more detailed analysis 
than simplex models. 

One additional point should be also discussed here. We have considered the model 
that is truly longitudinal (i.e., measurements at observation moments are correlated) but only 
examined and modeled data that are cross-sectional (i.e., different and uncorrelated samples 
are used for different moments). We solved the problem by considering only frequencies 
within state intervals at each time. However, the more general case might be studied for 
sufficiently large samples, viz. : when the data of interest were a matrix of transition 
frequencies (cell ij would contain the observed number of children changing from IQ state 
interval 1 at observed moment 1 to IQ state interval j at moment 2). The developed model 
would apply very nicely to this kind of data. 
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